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Abstract 

This  paper  focuses  on  the  development  of  parameter  estimation  techniques  for  models  quantifying  hysteresis 
and  constitutive  nonlinearities  in  ferroelectric  materials.  These  models  are  formulated  as  integral  equations  with 
known  kernels  and  unknown  densities  to  be  identified  through  least  squares  fit  to  data.  Due  to  the  compactness  of 
the  integral  operators,  the  resulting  discretized  models  inherit  ill-posedness  which  often  must  be  accommodated 
through  regularization.  The  accuracy  of  regularized  finite-dimensional  models  is  illustrated  through  comparison 
with  experimental  data. 
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1.  Introduction 

Piezoceramic  (PZT),  magnetostrictive  and  shape  memory  alloy  (SMA)  compounds  exhibit  hysteresis  and 
constitutive  nonlinearities  which  must  be  incorporated  in  models  and  model-based  control  algorithms  to  achieve 
the  novel  performance  capabilities  offered  by  the  materials.  A  number  of  modeling  strategies  for  these  compounds 
have  been  proposed  but  three  stand  out  in  the  sense  that  they  provide  unified  frameworks  for  characterizing 
hysteresis  in  ferroelectric,  ferromagnetic  and  ferroelastic  materials,  which  are  collectively  referred  to  as  ferroic 
compounds.  These  three  approaches  are  the  following:  (i)  homogenized  free  energy  models  [9,  13,  18],  (ii)  Preisach 
formulations  [1,  2,  12,  19],  and  (iii)  domain  wall  models  [6,  8,  11,  14,  15].  The  first  two  are  formulated  as  integral 
equations  whereas  the  domain  wall  models  are  typically  posed  as  differential  equations. 

To  simplify  the  discussion,  we  focus  on  the  estimation  of  parameters  in  the  homogenized  free  energy  model 
characterizing  the  hysteretic  field-polarization  relation  for  ferroelectric  materials.  It  is  demonstrated  in  [17]  that 
this  framework  can  be  used  to  characterize  hysteresis  in  general  ferroic  compounds  so  analogous  techniques  can 
be  employed  for  ferromagnetic  and  ferroelastic  materials.  It  is  illustrated  in  [16]  that  this  framework  provides  an 
energy  basis  for  certain  extended  Preisach  formulations,  and  compactness  results  analogous  to  those  established 
here  are  proven  for  classical  Preisach  operators  by  Iyer  and  Shirley  [7].  The  reader  is  also  referred  to  [5]  for 
details  regarding  the  recursive  identification  of  Preisach  density  functions. 

The  model  is  summarized  in  Section  2  and  compactness  of  the  integral  operator  is  established  in  Section  3. 
Parameter  identification  algorithms  are  summarized  in  Section  4  and  illustrated  in  Section  5  in  the  context  of 
characterizing  PZT5H  data. 
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2.  Polarization  Model 


It  is  demonstrated  in  [18]  that  the  hysteretic  relation  between  the  electric  field  E  and  polarization  P  in 
ferroelectric  materials  can  be  characterized  by  the  model 

pOC)  pOO 

[P{E)]{t)=  /  /  iy{E„Ee)[P{Ee  +  E;E„0]{t)dEedE,  (1) 

•Jo  •J  —  oo 

where  v  \s  &  material-dependent  parameter,  P  is  the  piecewise  linear  kernel  or  hysteron  depicted  in  Figure  1, 
and  Ec,Ee  respectively  denote  local  coercive  and  effective  fields.  The  kernel  has  the  general  form 

P{Ee  +  E;E„^)  =  -+PRS{E;E,,EO  (2) 

V 

where  ^  delineates  initial  dipole  orientations.  The  parameter  5  has  a  value  of  1  for  positively  oriented  dipoles 
and  -1  for  negative  orientations.  More  rigorously,  P  can  be  specified  as 

r  [P{E-,E,,Om  ,r(t)  =  0 

[P{E;Ec,^)](t)  =  I  ^-Pr  ,  r(t)  yf  0  and  F;(maxT(t))  = -Fie 

i  f  +  PR  ,  r(t)  yf  0  and  F;(maxT(t))  =  Flc 

where 

(  E-Pr  ,  F;(0)  <  -E, 

[P{E;  E,,  e)]  (0)  =  J  e  ,  -E,  <  E{0)  <  E, 

[  E+Pr  ,  E{0)  >  E, 

defines  the  initial  states  of  the  kernel  and 

r(t)  =  {te  {0,Tf]  I  E{t)  =  -E,  or  E{t)  =  E,] 
designates  the  set  of  switching  times. 

The  density  is  assumed  to  have  the  product  formulation  v{Ec,Ee)  =  vi{Ec)v2{E^)  where  vi  and  V2  satisfy 
the  physical  conditions 

(i)  vi{x)  defined  for  a;  >  0, 

(ii)  iy2{-x)  =  iy2{x), 

(3) 

(iii) 

|^^2(a:)|  <  C2e-“^" 

for  positive  ci, C2, 02.  The  goal  in  the  parameter  identification  problem  is  to  estimate  v  given  data  mea¬ 
surements  {{Ek,Pk)},k  = 


2 


By  invoking  the  physical  decay  criteria,  (1)  can  be  approximated  to  arbitrary  accuracy  by  consideration  of 


[P{E)]  (t)  =  JJ  E,)  [P{E  +  E,;E„  (t)  dE^  dE,  (4) 

0,2 

on  the  compact  domain 

^2  =  {{Ec,  Eg)  G  IR+  X  IR  I  v{Ec,Eg)  >  e}. 

Furthermore,  we  let  the  minimum  and  maximum  admissible  input  fields  be  denoted  Emin  and  Emax  and 
define 

—  \EminTEmax\’ 

We  consider  parameters  q  =  v  in  the  parameter  space 

Q  =  F"(f^2)  (5) 


and  define  the  observation  operator  CP  =  P{E)  on  the  observation  space 


P  (^Pmim  Pmax^  • 


(6) 


The  polarization  model  (4)  can  then  be  formulated  as 


y{E)  =  JCq{E) 


where 

E  e  C 

and  the  parameter-to-observation  operator  /C  is  defined  by 


1Cq  =  C 


k{-  +  Eg,  Ec)q{Eg,  Eg)  dEg  dEg. 


1^2 


(7) 


It  is  readily  observed  that  due  to  the  affine  construction  of  k  =  P,  k  G  L^{Ct)  and  k  G  L^(Cl)  where 

n  =  III  X  ^2  • 

The  property  that  k  G  is  typical  for  convolution  operators  whereas  k  G  facilitates  construction  of 

a  generalized  Fourier  basis  for  the  operator.  We  employ  this  latter  property  in  Section  3  to  establish  that  /C  is 
a  compact  operator. 

For  implementation  purposes,  it  is  demonstrated  in  [18]  that  Gaussian  quadrature  can  be  employed  to  ap¬ 
proximate  the  integrals,  thus  yielding  the  system 

Ni  Nj 

[PiE)]{t)  =  (8) 

i=i j=i 


where  Vi  and  Wj  denote  quadrature  weights  and  — !■  R.  To  formulate  (8)  as  a  linear  system,  we  define 

the  Ni  X  Nj  matrices  A{E)  and  4)  to  have  components 


[ME)hj  = 


E  +  Eg 
V 


+  PRS{E;Eg^,EgA 


V^W-j 


[Q]ij  =  v{Eg.,Eg-). 


For  N  =  Ni  ■  Nj,  we  define  the  fV  x  1  vector  q  and  1  x  fV  vector  a{E)  by 


q  =  vec((5)  ,  a{E)  =  [vec(xl(if))]^ 
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where  ‘vec’  denotes  the  vector  concatenation  of  the  respective  matrices.  The  discretized  polarization  model  (8) 
can  then  be  formulated  as  the  linear  system 


P{E)  =  a{E)q. 

We  note  that  77  is  considered  known  and  fixed  in  this  formulation  and  is  incorporated  in  a{E). 


(9) 


3.  Compactness  of  the  Polarization  Operator 

In  this  section,  we  establish  that  the  operator  JC  given  by  (7)  is  compact.  As  a  prelude,  we  state  the  following 
theorem  which  is  Theorem  5.24.8  from  [10]. 

Theorem  1.  Let  X  and  Y  he  Banach  spaces  and  let  K-n  '■  X  ^  Y,N  =  1,2, ,  he  a  sequence  of  compact  linear 
operators  converging  to  a  hounded  linear  operator  1C  :  X  ^  Y ;  that  is,  \\ICn  —  /C||  —>  0  as  ^  00.  Then  1C  is  a 
compact  linear  operator. 


Remark  1.  Consider  the  parameter  space  Q  and  observation  space  y  defined  in  (5)  and  (6).  The  integral 
operator  given  hy  (7)  is  then  a  compact  operator.  We  establish  this  hy  demonstrating  that  1C  is  the  limit  of  a 
sequence  of  finite  rank  operators  followed  by  the  use  of  Theorem  1. 


We  first  construct  an  orthonormal  basis  {fii}  for  Lf{Cl).  It  is  illustrated  in  [10]  that 


's/  ^max 


:  exp 


2^11  ■ 


s-  Er, 
l-E^ 


l=Q,±l,±2,--- 


forms  an  orthonormal  basis  for  L‘^{Cli). 
orthonormal  basis  for  Lfi{Cl)  is 


With  an  analogous  basis  definition  for  Lfi{Cl2),  it  follows  that  an 


film{s,t,v)  =  (pi{s)(pm{t,v) 


which  we  re-index  as  {4>i\. 

It  follows  that  every  /  G  Lf{Cl)  has  the  generalized  Fourier  series  representation 


i 


where  (•,  •)  denotes  the  usual  inner  product.  The  norm  representation 

ii/f  = 

i 

follows  from  Plancheral’s  theorem.  Moreover,  we  can  represent  K,  and  approximating  finite-rank  operators  /Cjv 

by 

ICf  =  Y,{f,fififii 

i 

N 

i=l 

where  fii  =  ICfii. 
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To  establish  the  convergence  /C  — >  /Cat,  we  note  that 

m-icNfw  = 


< 


< 


< 


E  I  a 1 11^*11 

z>A^+l 


E  E  11^* 


1/2 


1/2 


E  11^*1 


1/2 


where  the  third  inequality  follows  from  the  Schwartz  inequality.  Furthermore,  we  observe  that 


E  11^*11' 


E 


\ICUE)\^dE 

h 

E  J  jj  k{E  + Ee,Ec)(l)i{Ec,Ee)dEedE, 

^  02 

E  jJ  k{E  + Ee,Ec)4>i{Ec,Ee)dEedE, 

*  02 

JJ \k{E  +  Ee,E,)\‘^dEedE, 


dE 


dE 


dE  <  oo. 


The  last  step  follows  from  Plancheral’s  theorem.  The  convergence  of  llV'ilP  implies  that  X]i>Ar+i  ^  0 
as  iV  ^  oo.  Thus  for  £  >  0,  there  exists  such  that  for  N  >  N^, 


\\IC  -  /Cat II  =  sup 
f^o 


l|/C/ - /Cai/II 

ll/ll 


<  £ 


which  establishes  that 

lim  \\IC-ICn\\  =0. 

N—^oc) 

Since  the  range  of  /Cat  is  finite,  it  follows  that  /Cat  is  a  compact  operator.  The  compactness  of  1C  follows  from 
Theorem  1  since  it  is  the  norm  limit  of  a  sequence  of  compact  operators. 

The  compactness  of  K.  given  by  (7)  is  to  be  expected  since  it  is  a  special  case  of  a  Hilbert-Schmidt  operator 
which,  in  general,  can  be  characterized  as  having  an  kernel.  This  is  evidenced  by  the  fact  that  the  proof  given 
here  is  a  modification  of  that  in  [4]  for  Hilbert-Schmidt  operators  with  kernels  in 


4.  Parameter  Identification  Problem 


For  the  operator  K.  defined  in  (7),  data  P  corresponding  to  input  fields  E  G  {Emin,  E max),  and  parameter 

space  Q  =  L^{Cl2),  the  parameter  estimation  problem  can  be  formulated  as  follows:  find  q  G  Q  so  that 

ICq  =  P.  (10) 

We  note  that  (10)  has  a  classical  solution  if  and  only  if  P  G  TC{IC),  where  P(/C)  denotes  the  range  of  1C,  which, 
in  general,  will  not  be  true.  Instead  it  is  more  reasonable  to  consider  the  least  squares  problem 

mmT{q)  ,  T{q)  =  J\ICq-P\\y.  (11) 

q^Q  Z 
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However,  because  K,  is  compact  with  infinite  dimensional  range,  the  Moore-Penrose  inverse  K.^  is  discontinuous 
so  that  even  (11)  is  ill-posed  —  see  [3].  This  motivates  consideration  of  the  augmented  functional 

Tc.{q)  =  \\\ICq-Pfy  +  aJ{q)  (12) 

and  the  regularized  least  squares  minimization  problem 


minTaiq).  (13) 

qeQ 

The  regularization  parameter  a  >  0  controls  the  tradeoff  between  goodness  of  fit  to  the  data  and  stability 
whereas  the  penalty  functional  provides  stability  and  allows  the  inclusion  of  a  priori  information  regarding 
the  parameter  q.  One  choice  for  is  the  Tikhonov  functional  which  we  illustrate  in  the  context  of  the  discretized 
problem. 

To  formulate  the  finite-dimensional  parameter  estimation  problem,  we  modify  the  linearly  parameterized 
system  (9)  to  reflect  measured  data.  We  define  the  W  x  Nj  matrices 


Ek  +  Eg 
1 


+  Pn^iEk]  Ec^,  Ee-) 


ViWj 


[Q]ij  =  v{Ec,,Ee-) 


and  vector  concatenations 

q  =  vec(Q)  ,  Uk  =  [vec(Hfc)]^ 

so  that  q  and  Ofc  are  respectively  1  x  and  x  1  where  N  =  Ni  ■  Nj.  Additionally,  the  Nd  x  1  vectors  V  and 
V  are  defined  componentwise  by 

[V]k  =  P{Ek;q)  ,  [V]k  =  Pk  (14) 

and  the  Nd  x  N  matrix  A  is  defined  row-wise  by 


[A]k  =  ak- 

The  discretized  polarization  model  (8)  can  then  be  formulated  as  the  linearly  parameterized  system 


V{Ek)  =  Aq. 

The  unregularized  least  squares  problem  used  to  estimate  q  =  v  G  Q  =  given  measurements 

{{Ek,  Pk)},  k  =  1, Nd  is  the  following: 


mmT{q)  ,  T{q)  =  h\Aq-rf 

q£Q  Z 

subject  to  qi>t),  j  =  1, ...  ,N. 


(15) 


Here  ||  •  ||  denotes  the  Euclidean  norm  in  To  incorporate  Tikhonov  regularization,  we  consider  the  mini¬ 

mization  problem 


minT„(g)  ,  T{q)  =  i||Ag-'Pfcf  -k  ^||gf 

q£Q  Z  Z 


(16) 


subject  to  gi  >  0,  j  =  1,. . .  ,N. 


Techniques  for  choosing  a  to  avoid  oversmoothing  solutions  as  well  as  a  solution  algorithm  for  (16)  can  be  found 
in  Vogel  [20]. 
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5.  Validation  Example 

To  illustrate  attributes  of  the  least  squares  parameter  estimation  formulations  (15)  and  (16)  for  estimating 
the  N  =  Ni'  Nj  parameters  Eg^)},  we  consider  the  characterization  of  PZT5H.  The  unregularized  model 

fits  obtained  with  Ni  =  Nj  =  24  and  Ni  =  Nj  =  48  using  data  from  all  seven  hysteresis  loops  are  plotted 
in  Figure  2  whereas  those  obtained  using  the  same  quadrature  limits  in  the  regularized  functional  are  given  in 
Figure  3.  Without  regularization,  the  ill-posedness  associated  with  inversion  of  the  compact  integral  operator 
K,  yields  increasingly  inaccurate  model  predictions  as  discretization  limits  are  increased.  Regularization  through 
the  inclusion  of  the  penalty  term  stabilizes  the  pseudoinverse  by  shifting  singular  values  away  from  the 

origin  thus  yielding  the  highly  accurate  fit  observed  in  Figure  3. 


Figure  2:  PZT5H  data  and  model  fit  with  general  product  density  v  estimated  using  the  unregularized  functional 
(15)  with  data  from  all  7  loops,  (a)  Ni  =  Nj  =  24  {N  =  576),  and  (b)  Ni  =  Nj  =  48  {N  =  2304). 


Figure  3:  PZT5H  data  and  model  fit  with  general  product  density  v  estimated  using  the  regularized  Tikhonov 
functional  (16)  with  data  from  all  7  loops,  (a)  Ni  =  Nj  =  24  {N  =  576),  and  (b)  Ni  =  Nj  =  48  {N  =  2304). 


6.  Concluding  Remarks 

In  this  paper  we  have  developed  a  well-posed  framework  for  estimating  parameters  in  a  model  used  to 
characterize  the  hysteretic  and  nonlinear  field-polarization  relation  inherent  to  ferroelectric  materials.  The 
model  is  comprised  of  a  compact  integral  operator  with  infinite-dimensional  range  so  regularization  is  required 
to  guarantee  that  the  inverse  problem  is  well-posed.  It  is  illustrated  through  a  fit  to  PZT5H  data  that  Tikhonov 
regularization  yields  highly  accurate  model  fits  which  retain  stability  as  discretization  limits  are  increased. 
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